open(FH,"/Volumes/HD_2/Documents/lab/project/SNP/Normal_SNP.chr11.7M-9M.txt");
my @in=<FH>;
chomp(@in);
close(FH);
my $count = 0;
my $count_het = 0;
my $count_array_het = 0;
my $count_hom = 0;
my $fp = 0;
my $fn = 0;
my $total_count=0;
my $ag_count=0;
my $cg_count = 0;
my $count_array_hom = 0;
#open(FH,"/Users/yaping/Documents/workspace/gatk_read/Sting/snps.raw.variantonly.vcf");
open(FH,"/Volumes/HD_2/Documents/workspace/snps.raw.performance.newGATK.vcf");

my @vcf=<FH>;
chomp(@vcf);
close(FH);


my %vcf_hash = ();
foreach my $line(@vcf){
	next if $line=~/^\#/;
	my @splitvcf=split "\t",$line;
	my $key = $splitvcf[0]."~".$splitvcf[1];
	#next if $splitvcf[5] <=10;
	$vcf_hash{$key} = $line;
}
foreach my $line(@in){
	my @splitin=split "\t",$line;
	if($splitin[3] ne $splitin[4]){
		$count_array_het++;
		
	}
	else{
		$count_array_hom++;
	}
	my $key = "chr".$splitin[1]."~".$splitin[2];
	if( exists $vcf_hash{$key} ){
		$count++;
		my @splitvcf=split "\t",$vcf_hash{$key};
		if($splitvcf[9]=~/0\/1/ and $splitvcf[6]=~/PASS/){
			$count_het++;
			if($splitin[3] eq $splitin[4]){
				$fp++;
				print("fp:\t$splitin[0]\t$splitin[1]\t$splitin[2]\n");
			}
		}
		else{
			
				if($splitin[3] ne $splitin[4]){
					$fn++;
					if($line=~/\d+\tA\tB\tA\tG/ or $line=~/\d+\tA\tB\tG\tA/ or $line=~/\d+\tA\tB\tT\tC/ or $line=~/\d+\tA\tB\tC\tT/){
						$ag_count++;
					}
					print("false neg:\t$splitin[0]\t$splitin[1]\t$splitin[2]\n");
					
				}
				$count_hom++;
			
		}
		if($splitvcf[7]=~/CG/){
			print("CG site:\t$splitin[0]\t$splitin[1]\t$splitin[2]\n");
			$cg_count++;
		}
	}
	else{
		if($splitin[3] ne $splitin[4]){
			$fn++;
			if($line=~/\d+\tA\tB\tA\tG/ or $line=~/\d+\tA\tB\tG\tA/ or $line=~/\d+\tA\tB\tT\tC/ or $line=~/\d+\tA\tB\tC\tT/){
				$ag_count++;
			}
		
			print("not found-false neg:\t$splitin[0]\t$splitin[1]\t$splitin[2]\n");
		}
		$count_hom++;
		##print("$splitin[2]\n");
	}
	$total_count++;
}
my $fp_rate = $fp/$count_array_hom;
my $fn_rate = $fn/$count_array_het;
print("count is: $count\n");
print("total count is: $total_count\n");
print("het count is: $count_het\n");
print("array het count is: $count_array_het\n");
print("array hom count is: $count_array_hom\n");
print("false pos num is: $fp\t fp rate is: $fp_rate\n");
print("false neg num is: $fn\t fn rate is: $fn_rate\n");
print("ag count in fn is: $ag_count\n");
print("cpg count is: $cg_count\n");
